You are here: R workshop >> mapping demo

GBIF data access using R: Mapping demo

This session focuses on working with environment layers, mapping, cropping and masking layers using the Raster R-package and other tools. See also the Biostats making maps with R training resource.

GBIF data for king protea (Protea cynaroides) from South Africa

require('rgbif') # rOpenSci r-package for GBIF data
library('leaflet') # package to plot the Leaflet map
sp_name <- "Protea cynaroides"; kingdom <- "Plantae" # king protea, taxonKey=5637947
key <- rgbif::name_backbone(name=sp_name, kingdom=kingdom)$speciesKey
sp <- rgbif::occ_search(taxonKey=key, hasCoordinate=TRUE, country="ZA", limit=100)
sp_ZA = sp$data ### extracting the data table from the data frame
sp_ZA_m <- sp$data[c("name", "catalogNumber", "decimalLongitude","decimalLatitude", "basisOfRecord", "year", "stateProvince", "taxonKey", "occurrenceID")] ## Subset columns (for more useful map pop-up)
xy_ZA <- sp$data[c("decimalLongitude","decimalLatitude")] ## Extract coordinates
sp_ZA_xy <- sp$data[c("species", "decimalLongitude","decimalLatitude")] ## Input format for Maxent
leaflet() %>% addTiles() %>% addMarkers(lng = sp_ZA$decimalLongitude, lat = sp_ZA$decimalLatitude)
Map GBIF data for Protea cynaroides from South Africa
Map GBIF data for Protea cynaroides from South Africa

GBIF data for liverleaf (blåveis:no) from Norway

library('rgbif') # rOpenSci r-package for GBIF data
#library('mapr') # rOpenSci r-package for mapping (occurrence data) ## ERROR mapr deprecated
library('leaflet') # package to plot the Leaflet map
sp_name <- "Hepatica nobilis"; kingdom <- "Plantae" # liverleaf (blaaveis:no), taxonKey=5371699
key <- rgbif::name_backbone(name=sp_name, kingdom=kingdom)$speciesKey
sp <- rgbif::occ_search(taxonKey=key, hasCoordinate=TRUE, country="NO", limit=100)
sp_NO <- sp$data
sp_NO_m <- sp$data[c("name", "catalogNumber", "decimalLongitude","decimalLatitude", "basisOfRecord", "year", "municipality", "taxonKey", "occurrenceID")] ## Subset columns (for more useful map pop-up)
xy_NO <- sp$data[c("decimalLongitude","decimalLatitude")] ## Extract only the coordinates
sp_NO_xy <- sp$data[c("species", "decimalLongitude","decimalLatitude")] ## Input format for Maxent
### rgbif::gbifmap(sp, region = "norway") ### defunct gbifmap deprectated March 2023
### mapr::map_leaflet(sp_m, "decimalLongitude", "decimalLatitude", size=2, color="blue") ### defunct mapr deprectated March 2023
leaflet() %>% addTiles() %>% addMarkers(lng = sp_NO$decimalLongitude, lat = sp_NO$decimalLatitude)
Map GBIF data for Hepatica nobilis from Norway
Map GBIF data for Hepatica nobilis from Norway

Extract coordinates suitable for e.g. Maxent

library('utils')
xy <- sp$data[c("decimalLongitude","decimalLatitude")] ## Extract only the coordinates
sp_xy <- sp$data[c("species", "decimalLongitude","decimalLatitude")] ## Input format for Maxent
#utils::write.table(sp_xy, file="./mapping/sp_xy.txt", sep="\t", row.names=FALSE, qmethod="double")
#readLines("./mapping/sp_xy.txt", n=10)

Write example species occurrence data from R to file

#head(sp_ZA, n=5) ## preview first 5 records
utils::write.table(sp_NO, file="./mapping/sp_NO.txt", sep="\t", row.names=FALSE, qmethod="double") ## write to file
utils::write.table(sp_ZA, file="./mapping/sp_ZA.txt", sep="\t", row.names=FALSE, qmethod="double") ## write to file
## Extract coordinates
xy_NO <- sp_NO[c("decimalLongitude","decimalLatitude")] ## Extract coordinates
xy_ZA <- sp_ZA[c("decimalLongitude","decimalLatitude")] ## Extract coordinates

Read example species occurrence data into R

sp_ZA <- utils::read.delim("./mapping/sp_ZA.txt", header=TRUE, dec=".", stringsAsFactors=FALSE) ## read from file
xy_ZA <- sp_ZA[c("decimalLongitude","decimalLatitude")] ## Extract coordinates
sp_NO <- utils::read.delim("./mapping/sp_NO.txt", header=TRUE, dec=".", stringsAsFactors=FALSE) ## read from file
head(sp_ZA, n=5) ## preview first 5 records

Administrative borders for South Africa from GADM

#install.packages("geodata")
library(geodata)
gadm_ZA_2 <- gadm('ZAF', level=2, path="./data/gadm", resolution=2) ## level 0,1,2,...
plot(gadm_ZA_2, main="Adm. Boundaries South Africa Level 1")
points(xy_ZA, col='red', pch=20) ## plot species occurrence points to the map (smaller dots)
legend("top", title = "Legend", legend = "Protea cynaroides", pch = 20, pt.bg = "red", bty = "n")
GADM Norway
GADM Norway

Get administrative borders for Norway from GADM

Slow to render, many details for the coast…

#install.packages("geodata")
library(geodata)
gadm_NO_1 <- gadm('NOR', level=1, path="./data/gadm/", resolution=2) ## level 0,1,2,...
plot(gadm_NO_1, main="Adm. Boundaries Norway Level 1")
points(xy_NO, col='blue', pch=20) ## plot species occurrence points to the map (smaller dots)
legend("bottom", title = "Legend", legend = "H. nobilis", pch = 20, pt.bg = "blue", bty = "n")
GADM Norway
GADM Norway

Get world map from rworldmap

#install.packages("rworldmap")
library(rworldmap)
worldmap <- getMap(resolution = "coarse") ## getting a simple world map
plot(worldmap, col = "lightgrey", 
     border = "darkgray",
     xlim = c(-180, 180), ylim = c(-90, 90),
     bg = "aliceblue",
     asp = 1)
points(xy_NO, col='blue', pch=20, cex=1) # plot species occurrence points to the map
points(xy_ZA, col='red', pch=20, cex=1) # plot species occurrence points to the map
title("World map with some GBIF data points")

Border for Norway from maptool:wrld_simpl ## Extracting vector for one country (under development - TODO)

library(rworldmap)
map_NO <- gadm('NOR', level=0, path="./data/gadm/")
plot(map_NO)
points(xy_NO, col='blue', pch=20, cex=1) # plot species occurrence points to the map
#points(xy_ZA, col='red', pch=20, cex=1) # plot species occurrence points to the map
title("Country map with some GBIF data points")

Google API key – TESTING

library(ggmap)
## Since mid-2018, the Google Maps Platform requires a registered API key
## This is a PERSONAL and PRIVATE Google API key -- NB! this can generate costs
## Google requires a valid credit card to use its online cloud services
## So be sure to not share your API key -- do not save API keys in script files
##
## e.g. register_google(key = "mQkzTpiaLYjPqXQBotesgif3EfGL2dbrNVOrogg") (that's a fake key)
# this sets your google map for this session
# register_google(key = "[your key]")
##
## You can set it permanently by setting write = TRUE
## If you set it permanently it will be stored in your .Renviron file
## this sets your google map permanently:
# register_google(key = "[your key]", write = TRUE)
##

Base map from Google – Needs a Google API key

library(dismo)
library(ggmap)
##g_no <- gmap("Norway") ### Needs a Google API key
g_no <- ggmap("Norway") ### Needs a Google API key
trd <- geocode('Trondheim, Norway')
osl <- geocode('Oslo, Norway')
trd_merc <- Mercator(trd[, c('longitude', 'latitude')])
osl_merc <- Mercator(osl[, c('longitude', 'latitude')])
xy_merc <- Mercator(xy)
##
plot(g_no, interpolate=TRUE)
points(xy_merc, pch=20, col='blue') ## add species occurrence
points(trd_merc, pch='*', col='red', cex=3) ## add * for Trondheim
points(osl_merc, pch='*', col='red', cex=3) ## add * for Oslo
Base map from Google, with H. nobilis, Trondheim, Oslo
Base map from Google, with H. nobilis, Trondheim, Oslo

Testing remote sensing image for Trondheim downloaded from Landsat

I used the USGS LandsatLook Viewer and Sentinel2Look Viewer to download sattelite data for Trondheim.

## A rasterBrick is a multi-layer raster object (similar to a rasterStack)
rs_l <- brick('./data/landsat_trondheim_web_mercartor_wgs84.tif')
rs_s <- brick('./data/sentinel_trondheim_web_mercartor_wgs84.tif')
nlayers(rs_l); nlayers(rs_s) ## result = 3 layers
crs(rs_l); crs(rs_s) ## +proj=merc +a=6378137 +b=6378137 ... +units=m ...
ncell(rs_l); ncell(rs_s) ## rs_l = 100738 ## rs_s = 204530
dim(rs_l); dim(rs_s) ## rs_l = 241 418   3 ## rs_s = 362 565   3
res(rs_l);res(rs_s) ## rs_l = 30.09113 30.09113 ## rs_s = 20.03305 20.03305
#par(mfrow=c(2,1)) ## combining two plot with par(n_rows, n_columns) ## if landsat & sentinel
#plotRGB(rs_l, stretch="lin", axes=FALSE, main="Landsat True Color Composite") ## Landsat example
plotRGB(rs_s, stretch="lin", axes=FALSE, main="Sentinel True Color Composite") ## Sentinel example
Remote sensing data, sentinel, Trondheim
Remote sensing data, sentinel, Trondheim

GBIF data for taxon liverleaf (blåveis:no) from Trondheim

library('rgbif') # rOpenSci r-package for GBIF data
library('leafpop')
sp_name <- "Hepatica nobilis"; kingdom <- "Plantae" # liverleaf (blaaveis:no), taxonKey=5371699
key <- rgbif::name_backbone(name=sp_name, kingdom=kingdom)$speciesKey
bb <- c(10.2,63.3,10.6,63.5) # Trondheim
#bb <- c(5.25, 60.3, 5.4, 60.4) # Bergen
#bb <- c(18.7, 69.6, 19.2, 69.8) # Tromsoe
#bb <- c(10.6, 59.9, 10.9, 60.0) # Oslo
sp_bb <- rgbif::occ_search(taxonKey=key, hasCoordinate=TRUE, country="NO", geometry=bb, limit=100)
sp_bb_m <- sp_bb$data[c("name", "catalogNumber", "decimalLongitude","decimalLatitude", "basisOfRecord", "year", "municipality", "taxonKey", "occurrenceID")] ## Subset columns
leaflet() %>% addTiles() %>% addMarkers(lng = sp_bb$data$decimalLongitude, lat = sp_bb$data$decimalLatitude, popup = popupTable(sp_bb_m))
Map GBIF data with bounding box for Trondheim
Map GBIF data with bounding box for Trondheim

Further reading: